﻿using LxBasic;
using static LxLESS.MGTemp;

namespace LxCore
{
    public class AnisotropyDiffusionBase : EquationWithTransTerm
    {
        public MatrixR3C3 TransformMatrix = new DiagonalMatrixR3C3(1, 1, 1).ToMatrixR3C3();
        public override bool IsSymmetric()
        {
            return true;
        }
        public void SetZone_TransformMatrix(Vector3d axis1, Vector3d axis2)
        {
            if (axis1.Length() <= 0 || axis2.Length() <= 0)
            {
                LxLog.Error("旋转坐标不能为空");
                return;
            }
            axis1.Normalize();
            axis2.Normalize();
            if (Math.Abs(Vector3d.DotProduct(axis1, axis2)) > 1e-10)
            {
                LxLog.Error("坐标不垂直");
                return;
            }
            Vector3d axis3 = Vector3d.CrossProduct(axis1, axis2);
            axis3.Normalize();//可能不需要，但还是计算了
            TransformMatrix = MatrixR3C3.CalTransformMatrix(axis1, axis2, axis3);
        }

        public override void FormABOnFaces(ThreadInfo threadInfo)
        {
            {
                (int i1, int endp1) = Zone.Inner.FaceIdx.Split(threadInfo);
                var calpf = getPhysicalFaceCalc(Zone.Inner);
                for (; i1 < endp1; i1++)
                {
                    ref FvInnFace face = ref InnFaces[i1];
                    var pf = calpf.GetPhysicalFace(i1, ref face);
                    AInnO[i1] = -pf.gDiff;
                    InnFaceShare[i1] = pf.CalNonOrthogonalityCorrectionOfInner(ref face, Grad, InnX);
                }
            }
            for (int i1 = 0; i1 < EquationBounds.Length; i1++)
            {
                switch (EquationBounds[i1])
                {
                    case FluxWallBound fluxbou:
                        formAb_Bound_FluxCombined(Zone.Bounds[i1], fluxbou, threadInfo);
                        break;
                    case FixValueWallBound fix:
                        formAb_Bound_fixValue(Zone.Bounds[i1], fix, threadInfo);
                        break;
                }
            }
        }
        protected void formAb_Bound_fixValue(FvBound bound, FixValueWallBound fixBound, ThreadInfo threadInfo)
        {
            (int i1, int endp1) = bound.FaceIdx.Split(threadInfo);
            var calpf = getPhysicalFaceCalc(bound);
            var fixvalue = fixBound.FixValue.GetCellVarArr(bound);
            for (; i1 < endp1; i1++)
            {
                ref var face = ref BouFaces[i1];
                var pf = calpf.GetPhysicalFace(i1, ref face);
                double alun = ABouH[i1] = -pf.gDiff;
                double xbou = BouX[i1] = fixvalue.GetCellVar(i1);
                BouFaceShare[i1] = pf.CalNonOrthogonalityCorrectionOfBound(ref face, Grad, InnX, xbou) - alun * xbou;
            }
        }

        protected void formAb_Bound_FluxCombined(FvBound bound, FluxWallBound fluxBound, ThreadInfo threadInfo)
        {
            (int i1, int endp1) = bound.FaceIdx.Split(threadInfo);
            var calpf = getPhysicalFaceCalc(bound);
            int start = bound.FaceIdx.Start;
            if (fluxBound.IsConstSource())
            {
                var factors = fluxBound.BouFactors1;
                if (fluxBound.FixFlux_b.IsEqual0())
                {
                    Array.Clear(ABouH, i1, endp1 - i1);
                    Array.Clear(BouFaceShare, i1, endp1 - i1);
                    for (; i1 < endp1; i1++)
                    {
                        ref var face = ref BouFaces[i1];
                        var pf = calpf.GetPhysicalFace(i1, ref face);
                        double beforPhib = pf.gDiff;
                        double rightTerm = -pf.CalNonOrthogonalityCorrectionOfBound(ref face, Grad, InnX, BouX[i1]);
                        double B = rightTerm / beforPhib;
                        factors[i1 - start] = (B);
                    }
                }
                else
                {
                    Array.Clear(ABouH, i1, endp1 - i1);
                    var fb = fluxBound.FixFlux_b.GetFaceVarArr(bound);
                    for (; i1 < endp1; i1++)
                    {
                        ref var face = ref BouFaces[i1];
                        var pf = calpf.GetPhysicalFace(i1, ref face);
                        double flux_b = fb.GetFaceVar(i1);
                        flux_b *= face.Area;
                        double beforPhib = pf.gDiff;
                        double rightTerm = -pf.CalNonOrthogonalityCorrectionOfBound(ref face, Grad, InnX, BouX[i1]) + flux_b;
                        double B = rightTerm / beforPhib;
                        factors[i1 - start] = (B);//phi_b=A*phi_C+B
                        BouFaceShare[i1] = flux_b;
                    }
                }
            }
            else
            {
                var factors = fluxBound.BouFactors2;
                var fb = fluxBound.FixFlux_b.GetFaceVarArr(bound);
                var fa = fluxBound.FixFlux_a.GetFaceVarArr(bound);
                var fh = fluxBound.Fixh.GetFaceVarArr(bound);
                var ft = fluxBound.FixT.GetFaceVarArr(bound);
                for (; i1 < endp1; i1++)
                {
                    ref var face = ref BouFaces[i1];
                    var pf = calpf.GetPhysicalFace(i1, ref face);
                    double h = fh.GetFaceVar(i1);
                    double flux_a = fa.GetFaceVar(i1) - h;
                    double flux_b = fb.GetFaceVar(i1) + h * ft.GetFaceVar(i1);
                    flux_a *= face.Area;
                    flux_b *= face.Area;
                    double beforPhib = pf.gDiff - flux_a;
                    double rightTerm = -pf.CalNonOrthogonalityCorrectionOfBound(ref face, Grad, InnX, BouX[i1]) + flux_b;
                    double A = pf.gDiff / beforPhib;
                    double B = rightTerm / beforPhib;
                    factors[i1 - start] = (A, B);//phi_b=A*phi_C+B
                    BouFaceShare[i1] = flux_b + flux_a * B;
                    ABouH[i1] = flux_a * A;
                }
            }
            //if (fluxBound.FixFlux_a.IsEqual0() && fluxBound.FixFlux_b.IsEqual0())
            //{
            //    if (fluxBound.Fixh.IsEqual0())
            //    {
            //        Array.Clear(ABouH, i1, endp1 - i1);
            //        Array.Clear(BouFaceShare, i1, endp1 - i1);
            //    }
            //    else
            //    {
            //        var fh = fluxBound.Fixh.GetFaceVars(bound);
            //        var ft = fluxBound.FixT.GetFaceVars(bound);
            //        for (; i1 < endp1; i1++)
            //        {
            //            double h = fh.GetFaceVar(i1);
            //            formAb_Bound_FluxCombined(i1, -h, h * ft.GetFaceVar(i1));
            //        }
            //    }
            //}
            //else
            //{
            //    var fb = fluxBound.FixFlux_b.GetFaceVars(bound);
            //    var fa = fluxBound.FixFlux_a.GetFaceVars(bound);
            //    if (fluxBound.Fixh.IsEqual0())
            //    {
            //        for (; i1 < endp1; i1++)
            //        {
            //            formAb_Bound_FluxCombined(i1, fa.GetFaceVar(i1), fb.GetFaceVar(i1));
            //        }
            //    }
            //    else
            //    {
            //        var fh = fluxBound.Fixh.GetFaceVars(bound);
            //        var ft = fluxBound.FixT.GetFaceVars(bound);
            //        for (; i1 < endp1; i1++)
            //        {
            //            double h = fh.GetFaceVar(i1);
            //            formAb_Bound_FluxCombined(i1, fa.GetFaceVar(i1) - h, fb.GetFaceVar(i1) + h * ft.GetFaceVar(i1));
            //        }
            //    }
            //}
            //void formAb_Bound_FluxCombined(int faceID, double flux_a, double flux_b)
            //{
            //    ref var face = ref BouFaces[faceID];
            //    var pf = calpf.GetPhysicalFace(faceID, ref face);
            //    Vector3d smartV = pf.Normal - pf.gT * face.OToN;
            //    flux_a *= face.Area;
            //    flux_b *= face.Area;
            //    double svv = Vector3d.DotProduct(face.Normal, smartV) / InnCells[face.OSideCell].Volumn;
            //    double A_bb = (svv + pf.gT) - flux_a;
            //    double B_b = flux_b + (-Vector3d.DotProduct(Grad[face.OSideCell], smartV) + svv * BouX[faceID]);
            //    double A_bc_1 = pf.gT / A_bb;
            //    double B_b_1 = B_b / A_bb;
            //    factors[faceID - start] = (pf.gT / A_bb, B_b / A_bb);
            //    BouFaceShare[faceID] = flux_b + flux_a * B_b_1;
            //    ABouH[faceID] = flux_a * A_bc_1;
            //}
        }

        public override IConnectFactor GetConnectWallFactor(FvBound bound)
        {
            return new WallFactor(this, bound);
        }
        public class WallFactor : FactorCalcBase, IConnectFactor
        {
            public WallFactor(AnisotropyDiffusionBase equation, FvBound bound)
            {
                SetBasic(equation, bound);
                var matrixBound = equation.EquationBounds[bound.ID] as CoupledWallBound;
                this.fb = matrixBound.FixFlux_b.GetFaceVarArr(bound);
                this.fa = matrixBound.FixFlux_a.GetFaceVarArr(bound);
                this.calpf = equation.getPhysicalFaceCalc(bound);

            }
            IPhysicalFaceCalc calpf;
            IFaceVarArr<double> fb;
            IFaceVarArr<double> fa;

            public (double a11, double a12, double b1, double a21, double a22, double b2) GetFactor(int faceID)
            {
                ref var face = ref BouCFaces[faceID];
                var pf = calpf.GetPhysicalFace(faceID, ref face);
                double gammadiff = -pf.gDiff;
                double fcor = pf.CalNonOrthogonalityCorrectionOfBound(ref face, Grad, InnX, BouX[faceID]);
                return (-gammadiff, gammadiff, fcor,
                 gammadiff, -gammadiff - fa.GetFaceVar(faceID) * face.Area, fb.GetFaceVar(faceID) * face.Area - fcor);
            }

            public void SetFactor(int faceId, double A11, double B1)
            {
                BouA[faceId] = -A11;
                BouFaceShare[faceId] = B1;
            }
        }
        protected virtual IPhysicalFaceCalc getPhysicalFaceCalc(IIndexHost indexHost)
        {
            throw new NotImplementedException();
        }
        protected interface IPhysicalFaceCalc
        {
            PhysicalFace GetPhysicalFace(int i1, ref FvInnFace face);
            PhysicalFace GetPhysicalFace(int i1, ref FvBouFace face);
        }
    }
    public class OrthotropicDiffusion : AnisotropyDiffusionBase
    {
        public IField<double> k11;
        public IField<double> k22;
        public IField<double> k33;

        public void SetZone_DiffTerm(IField<double> k11, IField<double> k22, IField<double> k33)
        {
            this.k11 = k11;
            this.k22 = k22;
            this.k33 = k33;
        }
        protected override IPhysicalFaceCalc getPhysicalFaceCalc(IIndexHost indexHost)
        {
            return new OrthotropicFace()
            {
                k11 = k11.GetFaceVarArr(indexHost),
                k22 = k22.GetFaceVarArr(indexHost),
                k33 = k33.GetFaceVarArr(indexHost),
                TransformMatrix = TransformMatrix,
            };
        }
        protected class OrthotropicFace : IPhysicalFaceCalc
        {
            public IFaceVarArr<double> k11;
            public IFaceVarArr<double> k22;
            public IFaceVarArr<double> k33;
            public MatrixR3C3 TransformMatrix;
            public PhysicalFace GetPhysicalFace(int i1, ref FvInnFace face)
            {
                var pf = new PhysicalFace();
                var m = new DiagonalMatrixR3C3(k11.GetFaceVar(i1), k22.GetFaceVar(i1), k33.GetFaceVar(i1));
                var m2 = MatrixR3C3.Product(TransformMatrix, m);
                pf.SetFace(m2, face);
                return pf;
            }

            public PhysicalFace GetPhysicalFace(int i1, ref FvBouFace face)
            {
                var pf = new PhysicalFace();
                var m = new DiagonalMatrixR3C3(k11.GetFaceVar(i1), k22.GetFaceVar(i1), k33.GetFaceVar(i1));
                var m2 = MatrixR3C3.Product(TransformMatrix, m);
                pf.SetFace(m2, face);
                return pf;
            }
        }
    }
    public class AnisotropyDiffusion : AnisotropyDiffusionBase
    {
        public IField<double> k11;
        public IField<double> k12;
        public IField<double> k13;
        public IField<double> k21;
        public IField<double> k22;
        public IField<double> k23;
        public IField<double> k31;
        public IField<double> k32;
        public IField<double> k33;

        public void SetZone_DiffTerm(
            IField<double> k11, IField<double> k12, IField<double> k13,
            IField<double> k21, IField<double> k22, IField<double> k23,
            IField<double> k31, IField<double> k32, IField<double> k33)
        {
            this.k11 = k11; this.k12 = k12; this.k13 = k13;
            this.k21 = k21; this.k22 = k22; this.k23 = k23;
            this.k31 = k31; this.k32 = k32; this.k33 = k33;
        }
        protected override IPhysicalFaceCalc getPhysicalFaceCalc(IIndexHost indexHost)
        {
            return new AnisotropyFace()
            {
                k11 = k11.GetFaceVarArr(indexHost),
                k12 = k12.GetFaceVarArr(indexHost),
                k13 = k13.GetFaceVarArr(indexHost),
                k21 = k21.GetFaceVarArr(indexHost),
                k22 = k22.GetFaceVarArr(indexHost),
                k23 = k23.GetFaceVarArr(indexHost),
                k31 = k31.GetFaceVarArr(indexHost),
                k32 = k32.GetFaceVarArr(indexHost),
                k33 = k33.GetFaceVarArr(indexHost),
                TransformMatrix = TransformMatrix,
            };
        }
        protected class AnisotropyFace : IPhysicalFaceCalc
        {
            public IFaceVarArr<double> k11;
            public IFaceVarArr<double> k12;
            public IFaceVarArr<double> k13;
            public IFaceVarArr<double> k21;
            public IFaceVarArr<double> k22;
            public IFaceVarArr<double> k23;
            public IFaceVarArr<double> k31;
            public IFaceVarArr<double> k32;
            public IFaceVarArr<double> k33;
            public MatrixR3C3 TransformMatrix;
            public PhysicalFace GetPhysicalFace(int i1, ref FvInnFace face)
            {
                var pf = new PhysicalFace();
                var m = new MatrixR3C3(
                     k11.GetFaceVar(i1), k12.GetFaceVar(i1), k13.GetFaceVar(i1),
                     k21.GetFaceVar(i1), k22.GetFaceVar(i1), k23.GetFaceVar(i1),
                     k31.GetFaceVar(i1), k32.GetFaceVar(i1), k33.GetFaceVar(i1));
                var m2 = MatrixR3C3.Product(TransformMatrix, m);
                pf.SetFace(m2, face);
                return pf;
            }

            public PhysicalFace GetPhysicalFace(int i1, ref FvBouFace face)
            {
                var pf = new PhysicalFace();
                var m = new MatrixR3C3(
                     k11.GetFaceVar(i1), k12.GetFaceVar(i1), k13.GetFaceVar(i1),
                     k21.GetFaceVar(i1), k22.GetFaceVar(i1), k23.GetFaceVar(i1),
                     k31.GetFaceVar(i1), k32.GetFaceVar(i1), k33.GetFaceVar(i1));
                var m2 = MatrixR3C3.Product(TransformMatrix, m);
                pf.SetFace(m2, face);
                return pf;
            }
        }
    }
}
//protected override void UpdateBoundVars_Grad0(FvBound bound, ThreadInfo threadInfo)
//{
//    (int i1, int endp1) = bound.FaceIdx.Split(threadInfo);
//    var calpf = getPhysicalFaceCalc(bound);
//    for (; i1 < endp1; i1++)
//    {
//        ref var face = ref BouFaces[i1];
//        var pf = calpf.GetPhysicalFace(i1, ref face);
//        Vector3d smartV = pf.Normal - pf.gT * face.OToN;
//        ref double bv = ref BouX[i1];
//        double svv = Vector3d.DotProduct(face.Normal, smartV) / InnCells[face.OSideCell].Volumn;
//        double Ab = svv + pf.gT;
//        double Bb = -Vector3d.DotProduct(Grad[face.OSideCell], smartV) + svv * bv;
//        bv = (Bb + pf.gT * InnX[face.OSideCell]) / Ab;
//    }
//}